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A one-dimensional quantum spin model with the competing two-spin and three-spin interactions is investi- 
gated in the context of a tensor network algorithm based on the infinite matrix product state representation. The 
algorithm is an adaptation of Vidal's infinite time-evolving block decimation algorithm to a translation-invariant 
one-dimensional lattice spin system involving three-spin interactions. The ground-state fidelity per lattice site 
is computed, and its bifurcation is unveiled, for a few selected values of the coupling constants. We succeed in 
identifying critical points and deriving local order parameters to characterize different phases in the conventional 
Ginzburg-Landau- Wilson paradigm. 
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Introduction. Quantum phase transitions (QPTs) [01 0], 
occurring at absolute zero temperature, are driven by quan- 
tum fluctuations in quantum many-body systems. In con- 
densed matter, a theoretical description of such a system usu- 
ally involves two-body interactions. However, recent progress 
makes it possible to realize multi-body interactions with cur- 
rent technology in experiments, which, as theoretical analyses 
unveiled, produce exotic quantum phases. Examples include 
a spin model with competing two-spin and three-spin interac- 
tions in a system of trapped ions [ 3] and of ultracold atoms in 
triangular lattices [4] and cold polar molecules ]5|] . Unfortu- 
nately, the presence of multi-body interactions in a spin model 
normally renders it not exactly solvable. This limits the appli- 
cability of conventional analytic approaches. Therefore, it is 
desirable to develop approximated techniques or fully numer- 
ical approaches to address such a novel many-body system. 

In the past few years, significant progress has been made to 
develop efficient numerical algorithms to simulate quantum 
many-body lattice systems in the context of the tensor net- 
work representations [6, 7], resulting from advances in our un- 
derstanding of quantum entanglement present in ground-state 
wave functions for quantum many-body lattice systems [8]. 
In Ref. [6], Vidal introduced the infinite time-evolving block 
decimation (iTEBD) algorithm to find an approximate ground 
state for an infinite-size one-dimensional (ID) lattice system 
in the infinite matrix product state (iMPS) representation. In 
Ref. |7j], an infinite projected entangled-pair state (iPEPS) 
is proposed for an infinite-size quantum system in two and 
higher spatial dimensions. A peculiar feature of these al- 
gorithms is their ability to compute the ground-state fidelity 
per lattice site |@1- lfl3ll . Fidelity describes the distance be- 
tween two given quantum states, thus enabling us to capture 
drastic changes in ground-state wave functions for quantum 
many-body lattice systems undergoing QPTs H- JH. In 
Refs. ISlllOt], it was argued that the ground-state fidelity per 
lattice site is able to characterize QPTs, regardless of what 
type of internal order is present in quantum many-body states 
underlying various quantum lattice systems in condensed mat- 
ter. However, an adaptation of the algorithms is necessary to 
make it possible to investigate a spin model with multi-body 



interactions. 

In this work, we first adapt the iTEBD algorithm, making 
it suitable to generate a ground-state wave function for a ID 
translation-invariant spin model with the competing two-body 
and three-body interactions in the context of the iMPS rep- 
resentation. Then the ground-state fidelity per lattice site is 
computed, and its bifurcations is unveiled for different val- 
ues of the truncation dimension, with a bifurcation point as a 
pseudo-critical point, which tends to the critical point when 
the truncation is removed. In addition, the bifurcation implies 
symmetry spontaneous breaking, characterized by a local or- 
der parameter, a central concept in the conventional Landau- 
Ginzburg- Wilson paradigm [23, The latter is read off 
from a reduced density matrix of a representative ground-state 
wave function from a given phase. 

Infinite matrix product state algorithm adapted to a spin 
model with three-body interactions. We first describe an 
iMPS algorithm adapted to a spin model with three-body in- 
teractions. For an infinite-size translation-invariant quantum 
many-body system with three-body interactions described by 
a Hamiltonian H = h!- l ' l+1, ' +2 \ we introduce the iMPS rep- 
resentation that is translation-invariant under a three-site shift. 
That means we need three three-index tensors T' Alr , T ] Blr , and 
T k clr and three diagonal matrices Aa, Ab, and Ac to store a wave 
function in the context of the iMPS representation, as shown 
in Fig. [TJi). Here, i, j, and k are physical indices, which run 
over a ^-dimensional local Hilbert space, I and r are the inner 
bond indices, /, r =l,...,x, with^ being the truncation dimen- 
sion. This amounts to the statement that the unit cell of the 
iMPS representation consists of three consecutive sites. The 
ground-state wave function is projected out by resorting to the 
imaginary time evolution. 

The updating procedure of the algorithm is as follows. Ap- 
ply a three-site imaginary time evolution gate U (i, i+ 1 , i +2) = 
exp(-/i [ ''' +I '' +2]Ar ) over a time slice, At <k 1, to the unit cell of 
the iMPS representation, and contract all the tensors, as shown 
in Fig. [TJi)- As such, we have a tensor Mabc (Fig. |T](ii)). Re- 
shape the tensor Mabc into a xd 2 X xd matrix and perform 
the singular value decomposition (S VD) to the matrix, we get 
Mab, A~c, and Mc, as shown in Fig.[T](iii). Contract the tensor 
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FIG. 1 : (color online) The updating steps of an infinite matrix prod- 
uct state (iMPS) algorithm adapted to a spin model with three-body 
interactions. Here, and A, t (77 = A, B, C) are three-index tensors 
and diagonal matrices, respectively, (i) Apply a three-site gate Ujj,k 
to the cell, (ii) Contract all the tensors in (i) and reshape the tensor 
M ABC into a matrix, (iii) Perform the singular value decomposition 
(SVD) to yield M AB , A c , and M c . (iii) Contract the tensor M AB and 
A c and reshape it into a matrix, (iv) Perform the SVD to yield M A , 
A B , and M B . (v) Restore the iMPS representation with the identity 
resolution A n A~ [ = id inserted, (vi) The diagonal matrices A B , Ac, 
and tensors r,4, T B , and Tc are updated, as indicated by the dash-line 
ovals in (iv) and (v). 



Mab and Ac and reshape it into a yd X xd matrix and perform 
the SVD, we have Ma, Ab, and Mb, as seen in Fig. [T](iv). 
Insert the identity resolution /UAT 1 = id and contract Mq and 
A A l to yield fc- In Fig.Q](v), the identity resolution A n A n l = id 
(77 = A and C) is inserted, thus we are able to update the di- 
agonal matrices As, Ac, and tensors Fa, Tb, and Tc- Once this 
is done, we shift one site, apply the gate U(i + 1, i + 2, i + 3), 
and repeat the updating procedure. Afterwards, shift one site, 
apply the gate U(i + 2,i + 3, i + 4) and repeat the updating 
procedure again. This completes our updating procedure for 
a time slice At. Iterating until the ground-state energy is con- 
verged, the system's ground-state wave function is yielded in 
the iMPS representation. 

Model. We consider a ID spin model with the competing 
two-spin and three-spin interactions [3], which is described by 
the Hamiltonian: 

H = h 2 < rV < rV ' : ofo-f "trf 3 -hj]af, (1) 

i i i 

where J2 and J3 are, respectively, the two-spin and three-spin 
coupling constants, h is an external field along the x direction, 
and crjj ] (a — x,z) are the spin 1/2 Pauli operators at the i- 
th site. If the coupling constants are varied independently, 
it undergoes QPTs, with the occurrence of different quantum 
phases: an antiferromagnetic (AF) phase, a ferrimagnetic (F) 
phase, and a disordered paramagnetic (P) regime where spins 
are aligned along the x direction [3]. Note that, if J3 = 0, 
the model Hamiltonian H becomes quantum Ising model in 
a transverse field, which is exactly solvable, with the critical 
point at I1/J2 — 1. A schematic phase diagram [3] is shown in 
Fig. [2] Here, we fix h-1 and focus on three cases: (i) 72=0.4, 
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FIG. 2: (color online) Schematic phase diagram for a ID spin model 
with the competing two-spin and three-spin interactions. For Ji > 
and J3 > 0, there are three phases: a disordered paramagnetic (P) 
regime, an antiferromagnetic (AF) order phase, and a ferrimagnetic 
(F) phase. Three phases coexist in the tricritical point T. In addition, 
the model becomes quantum Ising chain in a transverse field if J3 = 
0. We focus on three typical cases: (i) The horizontal double arrow 
crosses the phase boundary from P to F along Jilh=QA. (ii) The 
horizontal dotted double arrow crosses a critical point along J 2 /h=0. 
(iii) The vertical double arrow crosses a discontinuous QPT point 
along J 3 /h=\.5. 



73 varies from = 0.80 to J3 = 1.00 (horizontal solid double 
arrow in Fig. [2]l. (ii) 72=0, 7? varies from 73 = 0.90 to 73 
1 . 1 (horizontal dotted double arrow in Fig. [2j». (iii) 73=1.5, 
72 varies from 72 = -0.1 to 72 = 0.1 (vertical solid double 
arrow in Fig.|2]i. 

Bifurcation in the ground-state fidelity per lattice site 
and symmetry spontaneous breaking. For two ground-state 
wave functions \(p(A\)) and \<p(A2)), with A\ and A2 being 
two values of the control parameter A, the ground-state fi- 
delity F(A\,A2) = \{(p(A\)\<p(A2))\ asymptotically scales as 
F{A\,A2) ~ d(A\,A2) L , where L is the number of the lat- 
tice sites. The scaling parameter, d(A\,A2), first introduced in 
Ref. IPJ], characterizes how fast the ground-state fidelity goes 
to zero when the thermodynamic limit is approached. Physi- 
cally, d{A\ , Az) is the averaged ground-state fidelity per lattice 
site, 

\nd(A x ,A 2 ) = lim F ^ Al \ (2) 

As discussed in Refs. |0,[l3l, it is well defined in the thermo- 
dynamic limit, although F(A\, A2) becomes trivially zero. The 
ground-state fidelity per lattice site is able to locate a QPT in 
a quantum many-body lattice system, regardless of what type 
of internal order is present. 

The iMPS representation allows to efficiently compute the 
ground-state fidelity per lattice site ifioll . Now we turn to the 
simulation results for the three cases, as selected above, of the 
model Hamiltonian (Q]i in the context of the iMPS algorithm. 

(i) 72 = 0.4 and h = 1, and 73 is chosen as the control 
parameter. The ground-state fidelity per lattice site is plotted 
in Fig.[3ja), with the truncation dimension^ =8, 16, and 32, 
respectively. Note that, 1^(73)), with J' 3 = 1.00, is chosen as 
a reference state. We observe that there is a bifurcation point 
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FIG. 3: (color online) Main: The ground-state fidelity per lattice site, 
d(J 3 , J' 3 ), as a function of J 3 , at a fixed J 3 , for a ID spin model with 
the competing two-spin and three-spin interactions. Here, Ji = 0.4 
(a) and J 2 = (b), respectively, and 73 is chosen as the control pa- 
rameter. In addition, \(p(J' 3 ))t with J' 3 = 1.00 and J' 3 = 1.10, is chosen 
as a reference state, respectively, in (a) and (b). We observe that there 
is a bifurcation point in the ground-state fidelity per lattice site, indi- 
cating that spontaneous symmetry breaking occurs when the control 
parameter J3 varies. The bifurcation point is a pseudo critical point 
for a given value of the bond dimension. Inset: An extrapolation 
with respect to the bond dimension is performed for the pseudo crit- 
ical points for ^=4, 8, 16, 32, and 64, yielding the critical points 
h c ~ 0.9034 (a) and J 3 ~ 1.0001 (b), when^- -> 00. 
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FIG. 4: (color online) The ground-state fidelity per lattice site, 
d{Ji, J' 2 ), as a function of J2, at a fixed J' 2 , for a ID spin model with 
the competing two-spin and three-spin interactions. Here, J3 = 1.5, 
and J 2 is chosen as the control parameter. In addition, the reference 
state \<p(J' 2 )) is fixed at J' 2 = 0.10. Spontaneous symmetry breaking 
occurs with J 2 > 0. A jump occurs in the ground-state fidelity per 
lattice site, meaning that the phase transition is discontinuous. The 
pseudo critical points are shown for different values of the truncation 
dimension, without any significant shift from^ = 8 to^ = 32. 



in the ground-state fidelity per lattice site, arising from spon- 
taneous symmetry breaking of the Z3 translational invariance 
under three-site shifts, when the control parameter J3 varies 
from J3 = 0.80 to ^3 = 1.00. The bifurcation point is a 
pseudo critical point for a given value of the truncation dimen- 
sion. The extrapolation of the pseudo critical points for ^=4, 
8, 16, 32, and 64 yields the critical point J 3c = 0.9034, when 



X — > 00. Therefore, a continuous QPT occurs at J 3c ~ 0.9034. 

(ii) J 2 — and h = 1, and ^3 is chosen as the control pa- 
rameter. The ground-state fidelity per lattice site, d(J^,J'^), 
as a function of J3, at a fixed J' v is plotted in Fig. Ob), 
with the truncation dimension x - 8, 16, and 32, respec- 
tively. The reference state \<p(J' 3 )) is chosen at J' 3 = 1.10. 
We observe that there is a bifurcation point in the ground- 
state fidelity per lattice site. Actually, four degenerate ground 
states arise from spontaneous symmetry breaking of the group 
G — Z2 x Z2 x Z2/Z2, when the control parameter J3 varies 
from 73 = 0.90 to J3 = 1.10. Here, four Z2 groups are, re- 
spectively, generated from ofcrf, crfcr*;, o~fcr^ an( j erf erf crF, 
where A, B, and C label three sites in the unit cell, as shown 
in Fig. Q] Note that a pseudo critical point J$ x is identified 
as a bifurcation point. Therefore, a continuous QPT occurs at 
Ji ~ 1.0001. 

(iii) The ground-state fidelity per lattice site, d( J2, J' 2 ), as a 
function of J2, at a fixed J' v is plotted in Fig. |4]for J3 = 1.5, 
with J2 as the control parameter. In addition, the reference 
state \ip{J' 2 )) is fixed at J' 2 = 0.1. Spontaneous symmetry 
breaking of the translation group Z3 under three-site shifts oc- 
curs for J 2 > 0. There is a jump at J2 = in the ground-state 
fidelity per lattice site, meaning that the phase transition is 
discontinuous. The pseudo critical points are shown for differ- 
ent values of the truncation dimension, without any significant 
shift, when^- is increased from^- = 8 to^- = 32. 

Local order parameters. In the conventional Ginzburg- 
Landau- Wilson paradigm, a central concept is local order pa- 
rameters, whose nonzero values characterize symmetry bro- 
ken phases. As advocated in Ref. [25], any local order pa- 
rameter may be derived from a reduced density matrix from a 
representative ground-state wave function in a symmetry bro- 
ken phase, if its iMPS representation is known. Actually, this 
follows from the non-zero-entry structure of a reduced den- 
sity matrix: in a symmetric phase, any reduced density matrix 
respects all the symmetry of the system, but in a symmetry 
broken phase, there are extra nontrivial entries in a reduced 
density matrix due to the fact that a symmetry operation is 
lost. Following this line of reasoning, we are able to derive 
local order parameters for different symmetry broken phases. 

In the Z3 symmetry broken phase, three degenerate ground 
states arise. The non-zero-entry structure of the one-site re- 
duced density matrix yields a local order parameter Omi = 
((of) - 2<of > + (of »/4, which is shown in Fig. Ha). In 
fact, we may equally choose any of the following three oper- 
ators as a local order parameter: (-2(0^) + (cr?) + (cr!r))/4, 
((of) - 2<of > + <of »/4, and «cr^> + (erf) - 2<of »/4, where 
A, B, and C label three sites in the unit cell, as shown in Fig.Q] 

In the G symmetry broken phase, four degenerate ground 
states arise. A local order parameter may be constructed from 
an analysis of the non-zero-entry structure of the one-site re- 
duced density matrix. One may choose (erf), (<rf ), and (cr^) 
as a local order parameter, which are able to distinguish four 
degenerate ground states. Their magnitudes take the same 
value, which is shown in Fig.|5jb). 

Now we take advantage of the finite-entanglement scaling 
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FIG. 5: (color online) Main: The local order parameter Ou\ = 
«cr^> - 2(crf) + (of»/4 shown in (a) and M2 = (cr z ) shown in 
(b) for a ID spin model with the competing two-spin and three-spin 
interactions, where Ji = 0.4 and Jj = 0, respectively. The compu- 
tation is performed with the truncation dimension^ = 8, 16 and 32. 
The pseudo critical points located by the order parameter are consis- 
tent with those obtained by the ground-state fidelity per lattice site. 
We see that M i and M i tend to saturate, with increasing^. Inset: 
The extrapolation of the semi-infinite chain von Neumann entropy is 
performed, yielding the central charge c ~ 0.796, with respect to the 
bond dimension x = 4, 6, 8, 10, 12, 32 and 64, in (a) and c ~ 0.994, 
with respect to the bond dimension x = 8, 16, 18, 24, 32, 36 and 50, 
in (b), very close to the exact results c = 4/5 and 1, respectively. 



of the von Neumann entropy to determine the nature of the 
criticalities. The von Neumann entropy S for a semi-infinite 
chain scales as S ~ cK/6\ogx, where c is the central charge, 
and k is the finite-entanglement scaling exponent 112611 . Our 
best fitting results, as shown in the insets of Fig.[5](a) and (b), 
yields the central charge c ~ 4/5 and c ~ 1, with k ~ 1.8077 
and k ~ 1 .2996, respectively. This implies that the transition 
points 73 e ~ 0.9034 and J^ c ~ 1.0001 are, respectively, in 
the same universality class as the three-state Potts universality 
class and the four-state Potts universality class. 

Conclusions. We have investigated a quantum spin model 
with the competing two-spin and three-spin interactions by 
exploiting a tensor network algorithm based on the iMPS rep- 
resentation. The algorithm itself is an adaptation of Vidal's 
iTEBD algorithm to a translation-invariant ID lattice spin sys- 
tem with three-spin interactions. This enables us to efficiently 
compute the ground-state fidelity per lattice site, which in turn 
makes it possible to locate critical points by identifying bifur- 
cation points in the ground-state fidelity per lattice site. For a 
few selected values of the coupling constants, we succeeded in 
identifying critical points, and deriving local order parameters 
to characterize different phases in the conventional Ginzburg- 
Landau- Wilson paradigm. 
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